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^ ' Abstract 

Sr ■ A novel method is presented for implementation of the extended 

^^ ■ mean-field theory incorporating two-body collisions. At a given time, 

C^ i stochastic imaginary time propagation of occupied states are used to 

generate a convenient basis. The quantal collision terms, including 
Q> ' memory effects, is then computed by a backward mean-field propa- 

gation of these single-particle states. The method is illustrated in 
an exactly solvable model. Whereas the usual TDHF fails to repro- 
Q ■ duce the long time evolution, a good agreement is found between the 

^ . extended TDHF and the exact solution. 
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1 Introduction 

Mean-field approximation is often employed to study static and dynamical 
properties of many-body systems in different branches of physics including 
atomic physics, condensed matter and nuclear physics |[l|. The complex 
quantal many-body dynamics is reduced to an effective one-body problem, 
in which particles move under an effective self-consistent mean-field poten- 
tial without experiencing correlations. Such an approach, which is usually 
referred to as Hartree-Fock (HF) in static limit and Time-Dependent Hartree- 
Fock (TDHF) in dynamics, is best suited at low energies at which binary col- 
lisions have little effects on dynamics since scattering into unoccupied states 
is severely inhibited due to Pauli blocking. The HF and TDHF models have 
been very successful in describing many static and dynamical properties in 
nuclear physics at low energies p|, |^, |^, |^. In heavy-ion collisions around 
Fermi energy, dynamical evolution exhibits strong dissipation and fluctuation 
properties. Mean-field alone is inadequate to describe the collision dynamics 
at these energies, and it is necessary to improve the one-body description 
beyond the mean-field approximation by incorporating binary collisions due 
to short range correlations into the equation of motion. This model is usu- 
ally referred to as the extended TDHF (ETDHF)[^, |^, ^, ^J. A large amount 
of work has been done on the formal development of the ETDHF, however, 
due to the numerical complexity often linked with conceptual problems, the 
applications of the theory on realistic situations remain a difficult problem. 



and only a few approximate quantal calculations exists so far ||T0|, 0, |T2|, [13| . 
Most of the applications of this theory have been carried out in semi-classical 
approximation, known as the Boltzmann-Uehling-Uhlenbeck (BUU) model 
I^ . It has been very successful for understanding a variety of features as- 
sociated with nuclear collisions at intermediate energies, including collective 
flow and particle production. It is generally taken for granted that the test 
particle simulation of the BUU model provides a good approximation for the 
average dynamics for small quantal systems that we face with in heavy ion- 
collisions. The actual test of the semi-classical models should be made by 
comparing the test particle simulations with quantum transport calculations. 
Such a comparisons has been made in the mean-field approximation in a re- 



cent work ||20|, |2T| , and demonstrated that the expansion dynamics exhibits 
quantum effects which persists up to rather high temperature T ^ 5MeV. 
Furthermore, the quantal features may play even more important role during 
the disassembly phase of the reaction. Therefore, it is off great interest to 



develop quantal transport descriptions of heavy-ion dynamics. 

In this paper, we propose a practical method for obtaining numerical 
solutions of the ETDHF theory in its quantal version, and illustrate the 
method in an exactly solvable model|jl5[. We study a fermionic case but the 
method can be applied to treat bosonic or even mixed systems. 

2 Extended Mean-Field Theory 

The basic idea of one body approaches is to project the dynamics onto the 
single-particle density matrix p{t). In the ETDHF theory the evolution of 
p(t) is determined by a transport equation, 

th^-[h[p],p]=K[p]. (1) 

The left hand side describes the evolution under the effective mean-field 
Hamiltonian which may be related to the energy functional E[p] as h[p] = 
dE[p]/dp [0, and the right hand side represents a quantal binary collision 
term ||l^ . It is convenient to express the collision term in the natural repre- 
sentation |^a(^)) which diagonalizes the single-particle density matrix 

p{t) = T.\Mt))n,{t) {Mt)\. (2) 



Then, the matrix elements of the collision term are given by [12 
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{Mt)\K [p] \iJy{t)) = -'- (F,,,,(t) + F;, ,(t)) (3) 



with 

FxAt) = I fdt'Y^{\5'\v^2\aP)^l {a'PWbS) ^\, (4) 

Z Jto 

{P^\'P5S'Paa'Pl3l3' — Paa' Pf3l3' P-yyPSS')\f, 

where pwit) = Sxy — Pwif) and fi2 denotes an effective residual interaction. 
The matrix elements of the residual interaction are given as 

{aP\v,2\X5)^ \, = {^^{t')^p{t')\ vr2 \iix{t')^&{t'))^ (5) 

where (■)^ denotes the antisymmetrized matrix element, and the states ip\{t') 
are obtained by backward propagation of the natural states from time t to 



t' by the mean-field propagator, ipx{t')) = t/^(t,t') j-i/'Ali^)), where the mean- 
field propagator is given by 

U{t, t') = Texp {-'- J^' h[p {s)]ds^ (6) 

with T as the time-ordering operator. In the collision term (4), px\'{t') = 

(V^a(^') \pi'^')\ 4'xfi't')) denotes the elements of density matrix at time t' in the 
basis of the backward propagated states and it can be expressed as 

PAA'(t') = E {'4^xit')\Mt'))nait') {Mt'MAt')) (7) 

where Uait') denotes the occupation number of the natural state lipait'))- 

The collision term involves usually two characteristic times: the corre- 
lation time Tcor of the matrix elements of the residual interactions and the 



relaxation time r^^ei of the occupation numbers ||18|| . The expression of the 
collision term is valid when the correlation time is much smaller than the 
relaxation time, Tcor <^ Treu which is usually referred to as the weak coupling 
limit. The time integration over the past history is evaluated over a time 
interval t — to which should be sufficiently larger than the correlation time 
Tcor- In this case, since during the correlation time, the natural states deviate 
only slightly from the states propagated with the mean-field, we may approx- 
imate the overlap as (ipxif) l^x'it')) ^ 6x\', and as a result, the collision 
term takes a simple form 

Fx,yit) ~ I f'dt'Y,{X6\v^2\af3)J^ {af3\v,2\X'S)^\, (8) 

Z Jto 

{nynsnaUp - nanpnxns)]^, ■ 

Furthermore, in accordance with the weak-coupling limit this expression may 
be further simplified by neglecting the evolution of the occupation numbers 
during the correlation time by taking them as nx{t') ~ nx{t). 

3 Numerical method 

For the purpose of numerical iteration, it is convenient to transform the 
transport eq.(l) into an integral form. Given the density matrix at a time 
t in terms of the occupied states \ipx{t)) and associated occupation numbers 



ft+At 

/ dsU{t + At,s)K{s)U{s,t + At) 



nx{t), its evolution during a time interval At may be expressed according to 

p{t + At) = U{t + At,t)p{t)U{t,t + At) (9) 

i rt+At 

where the first term represents the pure mean-field evolution and the second 
term is the perturbation caused by the collision term during the time interval 
At. Single-particle states t/j^it + At)) obtained by propagating the natural 
states from time t to t+At according to (6), provide a useful basis to calculate 
density matrix over short time intervals. In this representation, the elements 
of the density matrix at time t + At may approximately be given by 

1 rt+At 

pxx'it + At) ^ nxit)5xy - -^ / ds[Fxx>{s) + Fl,^{s)] (10) 

where -Faa'(s) is computed according to Eq. (4) or its approximate form 
(8) assuming that over the time interval At the natural states can be ap- 
proximated by tjjx{s)\. At time t + At the new natural states \^/Jx{t + At)) 
and their occupation numbers nx{t + At) are determined by diagonalizing 
the occupation matrix pxx'{t + At). Then, the iteration is continued into 
the next time step in a similar manner. The time interval At of the nu- 
merical iteration can be taken larger than the typical numerical time-step of 
the mean-field evolution, but it should be smaller than the correlation time 
to insure at most one collision takes place during the interval. When the 
collision term in eq.(9) is a small perturbation, a full diagonalization of the 
density matrix is not needed, and the new states and the changes in the oc- 
cupation numbers can be calculated using perturbation theory. However, to 
insure numerical stability, it is always possible to include higher orders. The 
occupation numbers are determined by the diagonal elements of (11) which 
can be transformed into a generalized master equation, 

^^-nxit) = f dt' {nx (f) Wt (t, t') - nx (f) Wx (t, t')} (11) 



dt 



to 



where the gain W^ and loss Wx kernels are given by 
1 



W+(t, t') = i E ^e {{)^S\v,2\af3)^ \t {af3\vu\XS)^ |a n, (t') n^ (f) ns (t')(12) 



and 

W,-(t, t') = ^Yl R^ {{\5\v^2\al3)j, \t (a/3|t;i2|A(5)^ 14 n, (t') np (f) n, (t') .(13) 



According to perturbation theory, when the states are non-degenerate, the 
new natural states at time t + At may be expressed as, 

\Mt + ^t)) = \^^{t + At)) (14) 

1 I — \ nt+At 

-72 E -7T. jrMyit + At)) / ds {F,y{s) + F:,,{s)}. 

h ^^nx[t) -nx'[t) I Jt 

This expression emphasizes importance of the non-diagonal terms in the colli- 
sion kernel which insure a proper transformation of the occupied states. With 
the diagonalization procedure the collision term itself decides the structure of 
the important single-particle states to be populated. For the application pre- 
sented in this paper, we have used a direct diagonalisation technics instead 
of a perturbative approach to define the waves functions. 

The collision term Fxx'{t) involves, in addition to the occupied single- 
particle states I'ip'x'^^it)) which are known, a complete set of unoccupied 
states. Specification of the most relevant set of intermediate states in a 
reasonable manner has been the major difficulty in the implementation of 
the ETDHF theory. Some effort has been made to generate the unoccupied 
states by the TDHF evolution in the same manner as the occupied states |T2| . 



Since it is difficult to guess from the beginning of the reaction, the relevant 
unoccupied states, which are strongly coupled with occupied states during a 
binary collision, may not be generated by a TDHF evolution. Moreover, the 
unoccupied states generated in this manner may escape to continuum before 
having a chance to participate in a possible binary collision. 

In this article, we propose an algorithm to construct the relevant unoccu- 
pied states at each stage of the iteration directly from the occupied states by 
employing a stochastic imaginary time propagator method as follows. At the 
beginning of each time step At, a series of states is generated by repeatec 
application of the imaginary time propagator f/i" = exp — /? [h[p] + 5h^'^^ 



r(«) 



,(«-!)/ 



,(")^ 



■,(0), 



on the occupied states Uf' Vx {t)) = K (t)} with ^r(^)/ = lV^r"(^))- 
The series of states obtained by application of the propagator exp{—j3h[p]) 
generate a restricted subspace which has the same symmetry properties as 
the occupied states. On the other hand relevant unoccupied states may 
have different symmetry properties, and therefore, may lie outside of this 
subspace. In order to ensure symmetry breaking and to perform a faster 
sampling of the relevant configuration |T9|, a stochastic part 6h^^^ is added to 
the the mean- field Hamiltonian at each application of the propagator Ug . 



In the case of open systems, the continuum may also be treated by adding a 
density-dependent constraining field into 5h^^^ . These series of linearly inde- 
pendent states together with the occupied states are orthonormalized using 
a Schmidt procedure after each application of U^ . In this manner, it is 
possible to construct a suitable subspace of properly orthonormalized single- 
particle basis {\'4'\{t))} including both occupied \'4''\'^^{t)) and unoccupied 
states \ip^°'^{t)). The subspace of the unoccupied states may be truncated 
further by diagonalizing the mean-field Hamiltonian h[p{t)] in this subspace 
and by removing the states with energies larger than a maximum amount. 
The high frequency components of the occupied states have only a minor 
effect in the collision term, since they damp out very quickly in time and 
because of the time integration in eq. (4) and (6). The construction with the 
help of the imaginary time propagator diminishes the high frequency com- 
ponents, and at the same time, it generates a space of the unoccupied states 
in terms of linear combinations of predominantly the low frequency compo- 
nents of each occupied states. In practice, for an approximate treatment of 
the collision term, it may be sufficient to generate the relevant unoccupied 
states by a few applications of Up on the occupied states. The sensitivity of 
the results on the approximate treatment can be checked by enlarging the 
subspace of unoccupied states. The unoccupied and occupied states are then 
evolved backward in time using the mean-field propagator order to properly 
take into account the past history in the collision term. We should note 
that, the high frequency components of the basis may become important for 
increasing available energy, and the method may depend on the structure of 
the residual interaction. Therefore, the method should be tested carefully 
before applications to nuclear collision dynamics around Fermi energy. 



4 First Application 

We illustrate the method in a one dimensional model problem of two identical 
fermions coupled to a total spin projection ms = 0, i.e. one with spin up and 
one with spin down, which are moving in an external anharmonic potential 
and are interacting via a short range two-body force. The Hamiltonian of 
the system is given by 

^ = E (it^ + ^(^^ - ^o)' + l(^^ - ^o)') + ^12 (15) 



with k = -0.04 MeV/fm2, k' = 0.08 MeV/fm^ xq = 9.3 fm. In this 
schematic model, we avoid the problem of construction of a G-matrix by em- 
ploying a simple force without a hard core. This force mimics the effective 
interaction usually employed in mean-field approaches, and also simplifies 
the discussion on the role of the collision term. We should note that the 
same effective interaction should be used in both the collision term and in 
the mean-field part of the ETDHF equation. The two-body interaction is 
taken as 

[Xi — X2) -, . , 

V12 = vo exp[ — J (16) 

with vq = —4 MeV and a = 2 fm. We consider the system initially prepared 
in a constrained equilibrium state at a temperature T=5 MeV. The initial 
two-particle density matrix is, therefore, represented by 

D,,(t = o) = E m'-^^^^^in (17) 

where the initial two-particle states are determined by solving a constrained 
Schroedinger equation, [H — \Q) |^j) = Ei |\E'j) with Q as a one-body con- 
straining field, which is taken as Q = J2i-^{xi){xi ~ ^oY with A = —0.24 
MeV/fm^ for Xi > xq and A = —0.12 MeV/fm^ for Xi < Xq. We follow the 
exact Liouville von Neumann evolution of the two-particle density matrix 
Di2(^) by evolving each two-particle states |^E'j(t)) with the time-dependent 
Schroedinger equation. We compare the exact evolution of the single-particle 
density matrix pi{t) = tr2Di2(t) with the approximate evolutions obtained 
in the TDHF and the ETDHF descriptions, starting with the same initial 
conditions. In the simulations, the numerical time-step is taken as 0.5 fm/c 
whereas the collision term is evaluated in larger time intervals of At = 6 
fm/c. In the calculation, we employ a harmonic potential with a stochastic 
strength for Sh and two iterations of the imaginary time propagator appeared 
to be sufficient to generate the relevant unoccupied states. Figure 1 shows 
the time evolution of the local single-particle density p{x, t) as a function 
of the position x. Until about 100 fm/c, the TDHF calculation shown by 
dashed lines provides a good approximation for the exact evolution indicated 
by dotted lines, but is not able to reproduce it for larger times. On the other 
hand, the ETDHF results indicated by solid lines are in good agreement 
with the exact evolution even at large times. Figure 2 illustrates the evolu- 
tion of the expectation value of the center of mass coordinate < X{t) >, as a 



function of time. Also for this observable, the ETDHF result (solid line) fol- 
lows closely the exact evolution (dotted line), whereas the TDHF calculation 
(dashed line) deviates from the exact evolution at large times. Time evolu- 
tion of the occupation numbers of time- dependent single-particle states are 
plotted in figure 3. The ETDHF calculations shown by dashed lines are very 
close to the exact evolution indicated by open circles, which are obtained by 
a direct diagonalization of pi{t) = tr2Di2(t). In the pure TDHF approach, 
the occupation numbers remain constant and equal to their initial values. 

5 conclusion 

Even though the formal development of the ETDHF theory was available for 
some time, only a few applications in some simplified model problems have 
been carried out so far. The major problem of the numerical implementation 
originates from the difficulties for a realistic treatment of the collision term 
in a suitable representation. In order to overcome this difficulty, we propose 
a possible method in which the most important unoccupied states, those 
strongly coupled with the occupied states through the collision term, are cal- 
culated dynamically at each time step and then evolved backward in time 
using the mean-field propagator in order to treat the past history in the col- 
lision term. We illustrate the method in an exactly solvable one-dimensional 
system of two-particles, and find that the description of the one-body den- 
sity matrix in the ETDHF is in good agreement with the exact evolution. 
In realistic applications to nuclear collision dynamics, the method may be 
implemented using the existing TDHF codes since it requires only standards 
calculation of the residual interaction matrix elements and usual mean-field 
propagation in real and imaginary time. In these applications, the compu- 
tational effort is not small, and may require a factor of about 20-100 more 
computation time than the corresponding TDHF simulation, but becomes 
manageable with high computational power of present day computer tech- 
nology. In fact, a first step in this direction has been recently taken |^2|]. In 



this work, collective vibrations at finite temperature have been investigated 
in the small amplitude limit of the ETDHF theory, including the coUisional 
term. Such a calculation requires the computation of many matrix element 
as in the full ETDHF description. The work for simulations of the full the- 
ory for large amplitude collective motion is currently in progress. We should 
also note that, the ETDHF theory is relevant to not only nuclear dynani- 



ics. Therefore, the proposed method may provide a useful tool for describing 
dissipation and fluctuation phenomena in other quantal systems. 
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Figure 1: Evolution of the local density. The exact calculations, the TDHF 
and the ETDHF results are shown by circles, dashed lines and solid lines, 
respectively. 
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Figure 2: Expectation value of the center of mass as a function of time. The 
exact calculation, the TDHF and the ETDHF results are shown by circles, 
dashed line and solid line, respectively. 
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Figure 3: Occupation numbers as a function of time. The exact calculations, 
the TDHF and the ETDHF results are shown respectively by circles, dashed 
and solid lines. 
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